Analysis of Colon Cancer Survival Data: Treatment Efficacy and Subgroup Identification
Vismaya S D | MSc Biostatistics| Manipal Academy of Higher Education
September 2025
1. Introduction and Objective
This report presents a comprehensive survival analysis of the
colon dataset, focusing on patients who experienced a death
event. The primary objective is to evaluate the efficacy of three
different treatment regimens: Observation (placebo), Levamisole (Lev),
and Levamisole combined with 5-FU (Lev+5-FU).
The analysis moves beyond a simple overall comparison to identify key prognostic factors and determine if the treatment effect is more pronounced in specific patient subgroups, providing actionable insights for clinical strategy.
2. Data Preparation and Exploration
Data Loading and Filtering
First, we load the necessary libraries and the colon
dataset. We then filter the data to include only records corresponding
to death events (etype == 2).
# Load and filter the data
colon <- survival::colon
colondeath <- filter(as_tibble(colon), etype == 2)
nrow(colondeath)## [1] 929## [1] 1858Baseline Patient Characteristics
A summary of the patient demographics across the three treatment arms shows a well-balanced distribution in terms of sample size, mean age, and gender. This balance gives us confidence that any observed differences in outcome are likely due to the treatment itself, rather than baseline imbalances.
# Summarize baseline characteristics by treatment group
colondeath %>%
group_by(rx) %>%
summarise(
N = n(),
`Mean Age` = round(mean(age, na.rm = TRUE), 1),
`Male (%)` = round(mean(sex == 1) * 100, 1)
) %>%
knitr::kable(caption = "Table 1: Baseline Patient Demographics by Treatment Arm")| rx | N | Mean Age | Male (%) |
|---|---|---|---|
| Obs | 315 | 59.5 | 52.7 |
| Lev | 310 | 60.1 | 57.1 |
| Lev+5FU | 304 | 59.7 | 46.4 |
3. Overall Survival Analysis
To get an initial understanding of treatment efficacy, we generate Kaplan-Meier survival curves for the three treatment arms.
# Fit and plot the overall survival curves
fit_km <- survfit(Surv(time, status) ~ rx, data = colondeath)
ggsurvplot(fit_km,
data = colondeath,
pval = TRUE, conf.int = TRUE,
legend.labs = levels(colondeath$rx),
palette = c("red", "blue", "green"),
title = "Overall Survival by Treatment (Colon Cancer Trial)"
)
Figure 1: Kaplan-Meier survival curves for the three treatment arms. The y-axis represents the probability of survival, and the x-axis represents time in days.
Interpretation
The Kaplan-Meier plot (Figure 1) shows a clear separation between the curves, with the Lev+5-FU (Green line) arm demonstrating the highest survival probability over the follow-up period. The log-rank test confirms this visual finding with a highly significant p-value (p < 0.05), indicating that the differences in survival among the three groups are statistically significant and not due to random chance.
4. Multivariable Risk Factor Analysis
While the Kaplan-Meier analysis is useful, it doesn’t account for other factors that might influence survival. We built a Cox Proportional Hazards model to assess the effect of treatment while adjusting for age, sex, and the number of positive lymph nodes.
Check Proportional Hazards Assumption - This is a critical step to validate the Cox model. We use the cox.zph() function, which tests the null hypothesis that the proportional hazards assumption is valid.
# Fit Cox model on full dataset
fit_cox_full <- coxph(Surv(time, status) ~ rx + age + sex + nodes, data = colon)
# Perform the test
ph_test <- cox.zph(fit_cox_full)
# Print the statistical results
# A p-value > 0.05 indicates the assumption holds for that variable.
print(ph_test)## chisq df p
## rx 1.03553 2 0.596
## age 0.00226 1 0.962
## sex 3.59295 1 0.058
## nodes 0.09340 1 0.760
## GLOBAL 4.81021 5 0.439# Create a visual plot of the Schoenfeld residuals
# If the assumption holds, the lines should be roughly flat and horizontal.
ggcoxzph(ph_test)
Interpretation of Your Results
GLOBAL Test (p = 0.439): This is the most important result. The overall test for your entire model has a p-value of 0.439, which is much higher than 0.05. This indicates that the model as a whole does not violate the proportional hazards assumption.
Individual Variables: All the individual predictors (rx, age, sex, nodes) also have p-values greater than 0.05. This means the effect of each variable on survival is constant over time, as required by the model
# FIX: To ensure perfect consistency and avoid plotting errors, we create a
# dedicated data frame with only the variables needed for this model.
# We also remove any rows with missing values for these variables.
# Fit Cox model on full dataset
fit_cox_full <- coxph(Surv(time, status) ~ rx + age + sex + nodes, data = colon)
# Forest plot works now
ggforest(fit_cox_full, data = colon)
Figure 2: Forest plot of the multivariable Cox model. Hazard Ratios less than 1 indicate a protective effect (reduced risk of death), while ratios greater than 1 indicate an increased risk.
Interpretation
The forest plot (Figure 2) provides several critical insights:
Treatment Effect: The Lev+5FU treatment is highly effective, reducing the hazard of death by 37% compared to observation (Hazard Ratio = 0.63, p < 0.001). The Levamisole-only treatment did not show a statistically significant benefit.
Prognostic Factors: The number of positive lymph nodes is the strongest predictor of outcome. For each additional positive node, the risk of death increases by 9% (HR = 1.09, p < 0.001). Age and sex were not significant predictors in this model.
Key Subgroup Insight: The benefit of the Lev+5FU treatment is most pronounced in high-risk patients (those with more than 3 positive nodes). The treatment provides a substantial survival advantage for this group, while its effect is minimal in low-risk patients.
5. Subgroup Analysis: Does Treatment Effect Vary by Risk?
The most advanced question is whether the treatment is more effective in certain patient subgroups. We hypothesize that the benefit of Lev+5-FU might be greater in higher-risk patients (those with more positive nodes). We tested this using an interaction term in our Cox model and then visualized the result by stratifying patients into “Low-Risk” and “High-Risk” groups based on the median number of nodes.
# Find the median number of nodes to define our groups
median_nodes <- median(colondeath$nodes, na.rm = TRUE)
# Create a new column to group patients
colondeath <- colondeath %>%
mutate(node_group = ifelse(nodes <= median_nodes, "Low Nodes", "High Nodes"))
# Split the data and fit separate models
low_nodes_data <- filter(colondeath, node_group == "Low Nodes")
high_nodes_data <- filter(colondeath, node_group == "High Nodes")
fit_low <- survfit(Surv(time, status) ~ rx, data = low_nodes_data)
fit_high <- survfit(Surv(time, status) ~ rx, data = high_nodes_data)
# Create a plot for each group
plot_low <- ggsurvplot(fit_low, data = low_nodes_data, title = "Low-Risk Group (<= 3 Nodes)", pval = TRUE, legend.labs = levels(colondeath$rx), palette = c("red", "blue", "green"))
plot_high <- ggsurvplot(fit_high, data = high_nodes_data, title = "High-Risk Group (> 3 Nodes)", pval = TRUE, legend.labs = levels(colondeath$rx), palette = c("red", "blue", "green"))
# Arrange the two plots side-by-side
plot_list <- list(plot_low, plot_high)
arrange_ggsurvplots(plot_list, print = TRUE, ncol = 2, nrow = 1, common.legend = TRUE, legend = "bottom")
Figure 3: Stratified Kaplan-Meier analysis. The plot on the left shows survival for low-risk patients (<= 3 nodes), while the plot on the right shows survival for high-risk patients (> 3 nodes).
Interpretation
The stratified analysis (Figure 3) reveals a powerful insight:
- In the Low-Risk Group (left), the survival curves are clustered together, and the p-value (p=0.3) is not significant. This means the choice of treatment has a minimal impact on outcomes for patients with less advanced disease.
- In the High-Risk Group (right), there is a dramatic separation of the curves. The Lev+5-FU treatment provides a substantial and statistically significant survival benefit (p < 0.05) compared to the other arms.
This confirms our hypothesis and is a key strategic finding: the therapy’s efficacy is most pronounced in the high-risk patient population.
6. Interactive Exploration
Finally, we can convert our overall survival plot into an interactive
dashboard using plotly. This allows for a more dynamic
exploration of the data, where users can hover to see specific data
points and toggle treatment arms on or off.
# Convert the static KM plot to an interactive one
library(plotly)
# Create static KM plot
km_plot_static <- ggsurvplot(fit_km,
data = colondeath,
pval = TRUE,
conf.int = TRUE,
legend.labs = levels(colondeath$rx),
palette = c("red", "blue", "green"),
title = "Interactive Survival by Treatment"
)
# Convert the ggplot part to interactive
ggplotly(km_plot_static$plot)7. Conclusion
This analysis demonstrates a clear and statistically significant survival benefit for the Lev+5-FU treatment regimen in colon cancer patients. The multivariable Cox model confirms this effect is robust after adjusting for key prognostic factors, with the number of positive lymph nodes being the strongest predictor of outcome.
Most importantly, the subgroup analysis reveals that the treatment’s benefit is primarily driven by its strong efficacy in high-risk patients (those with >3 positive nodes). This actionable insight is critical for optimizing clinical trial design and identifying the patient populations most likely to respond to this therapy.
Session Information
The following details show the R version and package versions used to generate this report, ensuring reproducibility.
## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_India.utf8 LC_CTYPE=English_India.utf8
## [3] LC_MONETARY=English_India.utf8 LC_NUMERIC=C
## [5] LC_TIME=English_India.utf8
##
## time zone: Asia/Calcutta
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] plotly_4.11.0 ggeffects_2.3.1 survminer_0.5.0 ggpubr_0.6.0
## [5] ggplot2_3.5.1 dplyr_1.1.4 survival_3.8-3
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.50 bslib_0.8.0 htmlwidgets_1.6.4
## [5] insight_1.4.0 rstatix_0.7.2 lattice_0.22-6 vctrs_0.6.5
## [9] tools_4.4.1 crosstalk_1.2.1 generics_0.1.3 tibble_3.2.1
## [13] pkgconfig_2.0.3 Matrix_1.7-0 data.table_1.16.4 RColorBrewer_1.1-3
## [17] lifecycle_1.0.4 compiler_4.4.1 farver_2.1.2 stringr_1.5.1
## [21] carData_3.0-5 htmltools_0.5.8.1 sass_0.4.9 yaml_2.3.10
## [25] lazyeval_0.2.2 Formula_1.2-5 pillar_1.10.1 car_3.1-3
## [29] jquerylib_0.1.4 tidyr_1.3.1 cachem_1.1.0 abind_1.4-8
## [33] km.ci_0.5-6 tidyselect_1.2.1 digest_0.6.37 stringi_1.8.4
## [37] purrr_1.0.2 labeling_0.4.3 splines_4.4.1 cowplot_1.1.3
## [41] fastmap_1.2.0 grid_4.4.1 cli_3.6.3 magrittr_2.0.3
## [45] broom_1.0.7 withr_3.0.2 scales_1.4.0 backports_1.5.0
## [49] rmarkdown_2.29 httr_1.4.7 gridExtra_2.3 ggsignif_0.6.4
## [53] zoo_1.8-12 evaluate_1.0.3 knitr_1.49 KMsurv_0.1-5
## [57] viridisLite_0.4.2 survMisc_0.5.6 rlang_1.1.5 xtable_1.8-4
## [61] glue_1.8.0 rstudioapi_0.17.1 jsonlite_1.8.9 R6_2.5.1